#include "Element_1D.h"
#include <cmath>

int main (int argc, char *argv[])
{
    const double PI = std::atan(1.0) * 4.0;
    const int N = 6;
    
    Point_1 mesh_pnts[N];
    //Set the mesh points
    mesh_pnts[0].set_coord(0.0);
    mesh_pnts[1].set_coord(PI / 4.0);
    mesh_pnts[2].set_coord(PI / 3.0);
    mesh_pnts[3].set_coord(PI / 2.0);
    mesh_pnts[4].set_coord(PI * 2.0 / 3.0);
    mesh_pnts[5].set_coord(PI);
    //set the mesh index
    for (int i = 0; i < N; i++)
	mesh_pnts[i].set_index(i);
    for (int i = 0; i < N; i++)
	std::cout << mesh_pnts[i] << ", ";
    std::cout << "\b\b" << std::endl;
    Mesh_1 mesh(mesh_pnts, N);
    Quad quad;
    quad.load_data("interval.quad");
    quad.list();
    double A[N][N];
    double f[N];
    for (int i = 0; i < mesh.get_n_grids(); i++)
	std::cout << mesh.get_element(i);
    std::cout << std::endl;
    const Quad_Info& q_info = quad.get_quad_info(10);
    std::cout << "accuracy: " << q_info.get_accuracy() << std::endl;
    int n_pnts = q_info.get_n_q_points();
    std::cout << "no. of quad points: " << n_pnts << std::endl;
    std::cout << "points:" << std::endl;
    double cont = 0.0;
    for (int i = 0; i < n_pnts; i++)
    {
	const double *pAw = q_info.get_pAw(i);
	double x = quad.l2g(0, PI, pAw[0]);
	cont += x * x * x * x * x * x * x * x * x * pAw[1] * 2.0;
	std::cout << x << ", " << pAw[1] << std::endl;
    }
    double J = quad.Jacobian(0, PI);
    std::cout << "int: " << cont * J << std::endl; 
    return 0;
};

